Trabajo individual de FAEDE+AEDM
Master en Bioinformática, Universidad de Murcia
Fernando Lucas Ruiz (fernando.lucas@um.es)
2024-02-26
Introducción
Este ejercicio está centrado en la búsqueda de perfiles de expresión génica en un conjunto de datos de microarrays de GEO(NCBI) que consta de 124 observaciones divididas en dos grupos/clases:
control (59) y
autism (65).
Las muestras de individuos autistas y de control se recogieron de personas de Phoenix. El grupo de cada observación se encuentra en el factor Class, columna 201 del fichero pone_0187371_s002.csv.
Librerias utilizadas
library(readr)
library(psych)
library(ggplot2)
library(viridis)
library(hrbrthemes)
library(tidyverse)
library(nortest)
library(car)
library(outliers)
library(corrr)
library(corrplot)
library(REdaS)
library(FactoMineR)
library(nFactors)
library(dendextend)
library(gplots)
library(RColorBrewer)
library(factoextra)
library(cluster)
library(MASS)Directorio de trabajo y carga de datos
En primer lugar, disponemos el directorio de trabajo para incorporar los datos del csv a nuestro entorno de RStudio.
setwd("/Users/fernandolucasruiz/Library/CloudStorage/OneDrive-UNIVERSIDADDEMURCIA/Documentos/Fernando/Master Bioinformatica/Asignaturas/Bioestadistica/tareas/JuanaMari")Análisis estadistico y limpieza de datos
Extraemos los datos de las columnas 103 a la 113 uniendo la columna de clase a mis datos.
Vemos que hay una diferencia muy grande en la expresion de las variables. Algunas variables tienen una pequeña discrepancia entre media y mediana por lo que podrÃan no ajustarse a una normal. Como era de esperar hay una n de 124 en todas las variables.
boxplot(data_flr[,1:11], col = colorRampPalette(c(rgb(255, 136, 0, maxColorValue = 255),
rgb(0, 94, 255, maxColorValue = 255)))(11))
Escalo los datos porque algunos valores entre variables son muy altos y otros muy bajos. De esta manera puedo comparar entre clases dentro de las variables y entre variables.
Ordeno las varianzas de mayor a menor.
data_flr <- as.tibble(scale(data_flr, center = T))
data_flr$class <- factor(data$Class)
data_flr <- data.frame(data_flr[, order(apply(data_flr, 2, var), decreasing = T)])
attach(data_flr)
boxplot(data_flr[,1:11], col = colorRampPalette(c(rgb(255, 136, 0, maxColorValue = 255),
rgb(0, 94, 255, maxColorValue = 255)))(11))
En esta tabla observamos los datos estadisticos descriptivos totales sin distinción de clase. Vemos que algunas variables tienen un skewness grande. PodrÃa deberse a una diferncia de clase (Control y Autista)
A continuación extraigo los datos estadisticos por clase (Control y Autismo) para ver si el skewness es debido a la clase.
Se intuyen diferencias en las medias y medianas entre los dos grupos.
En algunas variables se intuye que no son paramétricas debido a una diferencia entre media y mediana, además de un skewness distinto a cero. El caso más significativo es la variable 207084_at en los individuos autistas, que tiene un skewness alto y una Kurtosis alta por lo que tiene una distribucion platicúrtica.
## $Autism
## vars n mean sd median trimmed mad min max range skew
## X212864_at 1 65 -0.30 0.94 -0.52 -0.34 0.82 -2.31 2.05 4.36 0.45
## X208901_s_at 2 65 -0.29 0.99 -0.29 -0.29 1.12 -2.41 1.98 4.39 0.05
## X228949_at 3 65 -0.30 0.90 -0.46 -0.38 0.93 -1.52 1.82 3.34 0.71
## X215372_x_at 4 65 0.29 1.06 0.36 0.23 1.09 -1.91 4.33 6.24 0.81
## X227618_at 5 65 -0.29 0.93 -0.41 -0.33 0.81 -2.47 2.74 5.21 0.50
## X217142_at 6 65 0.30 0.90 0.31 0.30 0.94 -2.08 2.09 4.18 -0.15
## X207084_at 7 65 0.29 1.06 0.10 0.22 0.78 -1.91 5.77 7.69 2.03
## X1564876_s_at 8 65 0.29 1.04 0.03 0.25 0.88 -2.06 2.37 4.43 0.40
## X225363_at 9 65 -0.30 0.80 -0.28 -0.30 0.86 -2.29 1.33 3.62 -0.04
## X1563327_a_at 10 65 0.30 0.97 0.19 0.23 0.98 -1.58 2.80 4.38 0.56
## X231541_s_at 11 65 0.30 1.04 0.10 0.20 0.91 -1.51 3.66 5.17 0.90
## kurtosis se
## X212864_at -0.33 0.12
## X208901_s_at -0.52 0.12
## X228949_at -0.36 0.11
## X215372_x_at 1.77 0.13
## X227618_at 0.66 0.12
## X217142_at -0.41 0.11
## X207084_at 9.17 0.13
## X1564876_s_at -0.62 0.13
## X225363_at -0.52 0.10
## X1563327_a_at -0.16 0.12
## X231541_s_at 0.60 0.13
##
## $Control
## vars n mean sd median trimmed mad min max range skew
## X212864_at 1 59 0.33 0.97 0.27 0.29 0.73 -1.52 3.31 4.83 0.48
## X208901_s_at 2 59 0.32 0.92 0.29 0.32 1.10 -1.58 2.23 3.81 0.05
## X228949_at 3 59 0.33 1.01 0.27 0.29 1.04 -1.61 3.00 4.62 0.35
## X215372_x_at 4 59 -0.32 0.83 -0.31 -0.36 0.76 -2.15 1.87 4.02 0.44
## X227618_at 5 59 0.32 0.98 0.27 0.28 0.76 -1.56 2.93 4.49 0.58
## X217142_at 6 59 -0.33 1.01 -0.34 -0.35 0.95 -2.50 2.32 4.81 0.26
## X207084_at 7 59 -0.32 0.82 -0.37 -0.37 0.81 -1.87 2.50 4.37 0.68
## X1564876_s_at 8 59 -0.32 0.85 -0.47 -0.38 0.74 -2.02 1.72 3.74 0.58
## X225363_at 9 59 0.33 1.09 0.37 0.30 0.96 -2.08 2.97 5.05 0.15
## X1563327_a_at 10 59 -0.33 0.93 -0.40 -0.34 0.97 -2.67 1.78 4.45 0.09
## X231541_s_at 11 59 -0.33 0.85 -0.34 -0.37 1.04 -1.81 2.00 3.81 0.39
## kurtosis se
## X212864_at 0.48 0.13
## X208901_s_at -0.91 0.12
## X228949_at -0.28 0.13
## X215372_x_at 0.16 0.11
## X227618_at 0.21 0.13
## X217142_at 0.12 0.13
## X207084_at 1.11 0.11
## X1564876_s_at -0.24 0.11
## X225363_at -0.29 0.14
## X1563327_a_at -0.34 0.12
## X231541_s_at -0.41 0.11Distribución
Antes de ver las distribuciones de densidad, vamos a observar si hay ouliers que estén interfiriendo en los datos. Por ello creo boxplots para verlos. Observamos que en varios de estas variables existen valores ouliers.
par(mfrow=c(2,3))
for (i in 1:11) boxplot(data_flr[,i] ~ data_flr$class, xlab="", ylab = "", main = colnames(data_flr[i]), col = c("#F08080", "#4DD5F8"))
Decido reemplazar los outliers por los valores de su propia mediana mediante los argumentos fill = T, median = T. Tomo esta decisión porque no son muchas las variables y los valores en cada una de ellas.
data_flr_filtered <- rm.outlier(data_flr[,1:11], fill = T, median = T)
data_flr_filtered$class <- data_flr$class
data_flr <- data_flr_filteredVemos que siguen apareciendo outliers en alguno de ellos. Tomo la decisión de seguir porque están muy cerca de los rangos intercuartÃlicos. La señalización de la mediana en algunos casos no se sitúa en el centro por lo que intuimos que la distribución de los valores en algunos casos es no paramétrica.
par(mfrow=c(2,3))
for (i in 1:11) boxplot(data_flr[,i] ~ data_flr$class, xlab="", ylab = "", main = colnames(data_flr[i]), col = c("#F08080", "#4DD5F8"))
Vemos la distribución de las variables según la clase. Intuimos que hay diferencia entre medias y medianas según la clase por que no se solapan entre ellas. En la mayorÃa de distribuciones, vemos que tiene bastante cola en ambas direcciones ocn una distribución platicúrtica.
data_flr_gather <- data_flr %>%
gather(key = "text", value = "value", -class)
p <- data_flr_gather %>%
mutate(text = fct_reorder(text, value)) %>%
ggplot( aes(x=value, color=class, fill=class)) +
geom_density(alpha=0.6)+
#geom_histogram(alpha=0.6, binwidth = 5) +
ylim(c(0, 1)) +
scale_fill_viridis(discrete=TRUE) +
scale_color_viridis(discrete=TRUE) +
theme_ipsum() +
theme(
legend.position="none",
panel.spacing = unit(0.1, "lines"),
strip.text.x = element_text(size = 8)
) +
xlab("") +
labs(color = "Class") +
theme(legend.position = "bottom") +
guides(color = guide_legend(title = "Grupo", frame = TRUE)) +
facet_wrap(~text)
p
Contrastes
Realizo un estudio para ver si las variables son paramétricas o no. Selecciono las variables que no son paramétricas para realizar el contraste de medias apropiado en cada variable. Realizamos el test lillie ya que son más de 50 valores en cada variable y clase. Diferencio los ensayos por Control y Autista. Las variables X1564876_s_at, X1563327_a_at y X231541_s_at devuelven un pvalor significativo con el test Lillie por lo que sus valores son no paramétricos. Estos datos los guardo en un vector para la comparación de medias o medianas siguiente.
lista_noparametricos <- c()
for (i in colnames(data_flr)){
if (i == "class") {
break
}
sp_control <- lillie.test(data_flr[[i]][data_flr$class == "Control"])
if (sp_control$p.value < 0.05){
print(paste("La variable", i , "en las muestras Control tiene datos no paramétricos"))
lista_noparametricos <- append(lista_noparametricos, i)
}
sp_autism<- lillie.test(data_flr[[i]][data_flr$class == "Autism"])
if (sp_autism$p.value < 0.05){
print(paste("La variable", i , "en las muestras Autism tiene datos no paramétricos"))
lista_noparametricos <- append(lista_noparametricos, i)
}
}## [1] "La variable X1564876_s_at en las muestras Autism tiene datos no paramétricos"
## [1] "La variable X1563327_a_at en las muestras Autism tiene datos no paramétricos"
## [1] "La variable X231541_s_at en las muestras Autism tiene datos no paramétricos"En los qqplots podemos ver que las variables anteriormente predichas se alejan de la linea por lo que son no paramétricas. En el caso del Control X227618_at vemos que hay muchos valores que se dispersan pero realizando el test indidualmente da un pvalor de 0.1892. Con el test shapiro tampoco da significancia estadÃstica.
par(mfrow=c(2,3))
for (i in 1:11)
qqPlot(data_flr[,i][data_flr$class=="Control"], ylab= "",
main = paste("Control\n", colnames(data_flr[i])))
par(mfrow=c(2,3))
for (i in 1:11)
qqPlot(data_flr[,i][data_flr$class=="Autism"], ylab= "",
main = paste("Autism\n", colnames(data_flr[i]) ))
Con este bucle pretendo realizar test de contraste de medias o medianas si se da el caso de datos no paramétricos. En el caso de as variables no paramétricas realizo el test de U Mann-Witney.
En el caso de los datos paramétricos, primero realizo contraste de varianzas para observar si hay las varianzas son iguales. En el caso de si el test de contraste de varianza lance un pvalor < 0.05, realizamos el t.test con var.equal = F. Si no, con TRUE.
Como vemos, todas las variables tiene una diferencia significativa entre los individuos Control y los individuos con Autismo.
for (i in colnames(data_flr)){
if (i == "class") {
break
}
var_test <- var.test(data_flr[[i]] ~ data_flr$class)
if (i %in% lista_noparametricos){
wilcox_test <- wilcox.test(data_flr[[i]] ~ data_flr$class)
if (wilcox_test$p.value < 0.05){
print(paste("U Mann-Witney test significativo de",
i, "con un pvalor de:", round(wilcox_test$p.value, 6)))
}
} else {
if (var_test$p.value < 0.05){
t_test <- t.test(data_flr[[i]] ~ data_flr$class, var.equal = F)
if (t_test$p.value < 0.05){
print(paste("T test significativo (varianzas desiguales) de",
i, "con un pvalor de:", round(t_test$p.value, 6)))
}
} else {
t_test <- t.test(data_flr[[i]] ~ data_flr$class, var.equal = T)
if (t_test$p.value < 0.05){
print(paste("T test significativo (varianzas iguales) de",
i, "con un pvalor de:", round(t_test$p.value, 6)))
}
}
}
}## [1] "T test significativo (varianzas iguales) de X212864_at con un pvalor de: 0.000804"
## [1] "T test significativo (varianzas iguales) de X208901_s_at con un pvalor de: 0.001"
## [1] "T test significativo (varianzas iguales) de X228949_at con un pvalor de: 0.000814"
## [1] "T test significativo (varianzas iguales) de X215372_x_at con un pvalor de: 0.000741"
## [1] "T test significativo (varianzas iguales) de X227618_at con un pvalor de: 0.000942"
## [1] "T test significativo (varianzas iguales) de X217142_at con un pvalor de: 0.000746"
## [1] "T test significativo (varianzas iguales) de X207084_at con un pvalor de: 0.000421"
## [1] "U Mann-Witney test significativo de X1564876_s_at con un pvalor de: 0.000877"
## [1] "T test significativo (varianzas iguales) de X225363_at con un pvalor de: 0.00058"
## [1] "U Mann-Witney test significativo de X1563327_a_at con un pvalor de: 0.001329"
## [1] "U Mann-Witney test significativo de X231541_s_at con un pvalor de: 0.00245"Análisis de componentes principales
En primer lugar realizamos un estudio de correlación para ver si las variables se correlacionan entre si. En primer lugra mostramos la tabla con las correlaciones entre variables. En segundo lugar un plot que muestra el nivel de correlacion en la parte superior y unas elipses que muestra el nivel de correlación (más o menos eliptica) que existe entre esas dos variables y si es negativa o positiva dicha correlación según la direccion de la elipse.
Por ultimo, vemos la relación positiva o negativa de las variables.
En resumen, vemos que en la mayorÃa de pares de variables hay correlación entre ellas.
Para contrastar si nuestros datos son adecuados para hacer el analisis de componentes principales hay que pasar el test de esferecidad de Bartlett y el test de Kaiser-Meyer-Olkin.
En ambos tests rechazamos la hipotesis nula por lo que nuestros datos son adecuados para realizar el ACP. Además el test KMO nos da un criterio de 0.77 que es un criterio alto de adecuación para el ACP.
## Bartlett's Test of Sphericity
##
## Call: bart_spher(x = data_flr[, 1:11])
##
## X2 = 382.914
## df = 55
## p-value < 2.22e-16##
## Kaiser-Meyer-Olkin Statistics
##
## Call: KMOS(x = data_flr[, 1:11])
##
## Measures of Sampling Adequacy (MSA):
## X212864_at X208901_s_at X228949_at X215372_x_at X227618_at
## 0.7831949 0.7998157 0.8279270 0.7545901 0.7556934
## X217142_at X207084_at X1564876_s_at X225363_at X1563327_a_at
## 0.8264883 0.7861366 0.7863138 0.7781170 0.7236957
## X231541_s_at
## 0.5954110
##
## KMO-Criterion: 0.7691112El siguiente paso es realizar el análisis per se de los componentes principales de nuestros datos. Para ello escogemos 11 componentes ya que tenemos 11 variables.
En los datos podemos estadisticos podemos ver que seguramente escojamos 3 componentes ya que el cuarto componente tiene menos de 1 varianza en su componente. No llega al 80 % de las variables pero tenemos que trabajar con los datos que mejor nos convengan. Entre las 3 componentes recogemos el 57.8 de las varianza, como vemos en los datos y la gráfica screeplot.
En estos datos podemos la contribución de individuos y variables en las sucesivas componentes principales. Más adelante entraremos en profundidad en cada uno de ellos.
##
## Call:
## PCA(X = data_flr, scale.unit = T, ncp = 11, quali.sup = 12, graph = F)
##
##
## Eigenvalues
## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5 Dim.6 Dim.7
## Variance 3.837 1.447 1.077 0.949 0.802 0.750 0.606
## % of var. 34.878 13.154 9.795 8.625 7.287 6.814 5.512
## Cumulative % of var. 34.878 48.032 57.827 66.452 73.739 80.553 86.064
## Dim.8 Dim.9 Dim.10 Dim.11
## Variance 0.481 0.454 0.341 0.256
## % of var. 4.375 4.129 3.103 2.329
## Cumulative % of var. 90.439 94.568 97.671 100.000
##
## Individuals (the 10 first)
## Dist Dim.1 ctr cos2 Dim.2 ctr cos2 Dim.3
## 1 | 3.905 | 2.992 1.881 0.587 | -1.303 0.947 0.111 | 0.903
## 2 | 2.753 | 1.562 0.513 0.322 | 0.200 0.022 0.005 | 0.044
## 3 | 1.962 | 1.362 0.390 0.482 | -0.321 0.057 0.027 | -0.004
## 4 | 2.464 | 0.698 0.103 0.080 | -2.093 2.441 0.721 | -0.260
## 5 | 1.911 | 1.089 0.249 0.325 | 0.105 0.006 0.003 | 0.548
## 6 | 4.186 | -0.326 0.022 0.006 | -2.425 3.278 0.336 | -0.724
## 7 | 1.796 | -0.006 0.000 0.000 | 0.113 0.007 0.004 | 0.899
## 8 | 3.333 | 1.729 0.628 0.269 | -0.538 0.161 0.026 | 0.481
## 9 | 1.649 | 0.883 0.164 0.287 | 0.707 0.278 0.184 | 0.810
## 10 | 3.489 | 0.042 0.000 0.000 | -1.735 1.677 0.247 | 0.400
## ctr cos2
## 1 0.610 0.053 |
## 2 0.001 0.000 |
## 3 0.000 0.000 |
## 4 0.050 0.011 |
## 5 0.225 0.082 |
## 6 0.392 0.030 |
## 7 0.605 0.251 |
## 8 0.173 0.021 |
## 9 0.491 0.241 |
## 10 0.120 0.013 |
##
## Variables (the 10 first)
## Dim.1 ctr cos2 Dim.2 ctr cos2 Dim.3 ctr
## X212864_at | 0.744 14.443 0.554 | 0.374 9.662 0.140 | -0.093 0.803
## X208901_s_at | 0.665 11.513 0.442 | -0.030 0.061 0.001 | 0.241 5.374
## X228949_at | 0.522 7.095 0.272 | 0.301 6.256 0.091 | 0.092 0.783
## X215372_x_at | -0.601 9.410 0.361 | 0.026 0.045 0.001 | 0.588 32.130
## X227618_at | 0.669 11.678 0.448 | 0.340 7.979 0.115 | 0.297 8.211
## X217142_at | -0.507 6.709 0.257 | 0.277 5.309 0.077 | -0.027 0.067
## X207084_at | -0.321 2.689 0.103 | 0.358 8.863 0.128 | 0.438 17.772
## X1564876_s_at | -0.700 12.771 0.490 | 0.239 3.934 0.057 | 0.355 11.719
## X225363_at | 0.738 14.209 0.545 | 0.383 10.163 0.147 | 0.138 1.772
## X1563327_a_at | -0.479 5.975 0.229 | 0.527 19.217 0.278 | -0.379 13.339
## cos2
## X212864_at 0.009 |
## X208901_s_at 0.058 |
## X228949_at 0.008 |
## X215372_x_at 0.346 |
## X227618_at 0.088 |
## X217142_at 0.001 |
## X207084_at 0.191 |
## X1564876_s_at 0.126 |
## X225363_at 0.019 |
## X1563327_a_at 0.144 |
##
## Supplementary categories
## Dist Dim.1 cos2 v.test Dim.2 cos2 v.test Dim.3
## Autism | 0.945 | -0.916 0.939 -5.442 | 0.171 0.033 1.652 | 0.006
## Control | 1.041 | 1.009 0.939 5.442 | -0.188 0.033 -1.652 | -0.007
## cos2 v.test
## Autism 0.000 0.069 |
## Control 0.000 -0.069 |fviz_screeplot(data_flr.pca, addlabels = T, ylim = c(0,40), barfill = "dodgerblue3", linecolor = "#F08080")
Estudio variables
En primer lugar vamos a estudiar el papel de las variables en las componentes principales.
En la gráfica inferior podemos ver las coordenadas de las variables en las componentes principales 1 y 2. Vemos que X225363_at y X212864_at tienen un alto grado de contribución para ambas componentes ya que se situa en los 45º con respecto a las dos dimensiones.
fviz_pca_var(data_flr.pca, col.var="contrib",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = TRUE
)
fviz_pca_var(data_flr.pca, col.var="contrib",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = TRUE,
axes= c(1,3)
)
En esta gráfica podemos ver el grado de contribución de cada variable por componente principal. Como he dicho anteriormente, estudio los 3 primeros componentes. Como hemos visto en el mapa anterior, se corrobora que X225363_at y X212864_at tienen un alto grado de contribución en las dos primeras componentes pero no en la tercera componente. En cambio, la tercera componente viene marcada por la contribucion de X215372_x_at.
Individuos
Como vemos en la gráfica, los individuos presentes en la zona donde se cortan los ejes tienen una contribución muy baja en estas componentes principales.
Por otro lado, tenemos un alto grado de contribución del individuo 80 (izquierda eje X) en la componente 1 o el individuo 46 en la componente 2.
En la gráfica de barras se muestra la contribución de cada individuo en cada componente. A modo de resumen, el individuo 80 es el que más contribuye a la componente 1, el individuo 21 a la componente 2 y el individuo 69 a la componente 3.
fviz_pca_ind(data_flr.pca, col.ind = "cos2",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = TRUE
)
fviz_pca_ind(data_flr.pca, col.ind = "cos2",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = TRUE,
axes = c(1,3)
)
Conjunto
En estas gráficas mostramos el overlapping de individuos y variables.
Podemos decir que hay mucho overlap entre la distribución de individuos controles y autistas, pero parece que algunas variables influyen más en los individuos autistas y otros más en los controles.
Por ejemplo X215372_x_at parece ser que tiene que ver con los individuos autistas ya que se situa a la izquierda de la componente 1. En cambio, X208901_s_at está más presente en controles.
fviz_pca_biplot(data_flr.pca, col.ind = class, habillage = data_flr$class, addEllipses = T, axes = c(1,3))
Elección de componentes principales
El resultado optimo del número óptimo de componentes principales es 2, pero voy a escoger 3 componentes principales porque considero que hay pocos datos para seguir hacia delante con el ensayos de clusterización. También lo haré sin reducción para ver si mejora la jerarquización.
dim_data_flr <- nrow(data_flr)
ev <- eigen(cor(data_flr[,1:11]))
ap <- parallel(subject=dim_data_flr,var=ncol(data_flr[,1:11]),rep=100,cent=.95)
nS <- nScree(ev$values, ap$eigen$qevpea)
vp <- min(nS$Components) + 1
print(paste("El número idoneo de componentes principales es", vp))## [1] "El número idoneo de componentes principales es 2"Análisis jerárquico
Vamos a realizar el ensayo para una jerarquización no supervisado de los datos.
En primer lugar hacemos la matriz de los datos compilados por componentes principales del ensayo anterior y de los datos totales.
Ahora realizamos la clusterización mediante 5 métodos de enlace como son completo, simple, mediante el promedio, centroide o mediante enlace Ward.
Podemos ver que el método de enlace con mejor puntuación en el caso de todas las variables es el que mide la clusterización mediante la media con 0.58.
hclust_data_flr_complete <- hclust(dist_data_flr, method = "complete")
hclust_data_flr_avg <- hclust(dist_data_flr, method = "average")
hclust_data_flr_single <- hclust(dist_data_flr, method = "single")
hclust_data_flr_centroid <- hclust(dist_data_flr, method = "centroid")
hclust_data_flr_ward <- hclust(dist_data_flr, method = "ward.D")
complete <- cor(dist_data_flr.comp, cophenetic(hclust_data_flr_complete))
average <- cor(dist_data_flr.comp, cophenetic(hclust_data_flr_avg))
single <- cor(dist_data_flr.comp, cophenetic(hclust_data_flr_single))
centroid <- cor(dist_data_flr.comp, cophenetic(hclust_data_flr_centroid))
ward <- cor(dist_data_flr.comp, cophenetic(hclust_data_flr_ward))
best <- data.frame(complete, single, centroid, average, ward)
round(best, 2)Podemos ver que el método de enlace con mejor puntuación en el caso de los variables compiladas es el que mide la clusterización mediante la media. La compilación mejora el grado de bondad del clusterización.
hclust_data_flr_complete.comp <- hclust(dist_data_flr.comp, method = "complete")
hclust_data_flr_avg.comp <- hclust(dist_data_flr.comp, method = "average")
hclust_data_flr_single.comp <- hclust(dist_data_flr.comp, method = "single")
hclust_data_flr_centroid.comp <- hclust(dist_data_flr.comp, method = "centroid")
hclust_data_flr_ward.comp <- hclust(dist_data_flr.comp, method = "ward.D")
complete <- cor(dist_data_flr.comp, cophenetic(hclust_data_flr_complete.comp))
average <- cor(dist_data_flr.comp, cophenetic(hclust_data_flr_avg.comp))
single <- cor(dist_data_flr.comp, cophenetic(hclust_data_flr_single.comp))
centroid <- cor(dist_data_flr.comp, cophenetic(hclust_data_flr_centroid.comp))
ward <- cor(dist_data_flr.comp, cophenetic(hclust_data_flr_ward.comp))
best <- data.frame(complete, single, centroid, average, ward)
round(best, 2)Buscando el número de clusters generados, podemos ver que tanto en los datos completos como en los datos compilados según componentes principales, vemos que se generan 3 clusters. En las tablas vemos que la discrepancias aparece en los clusters 2 y 3 pero en el 1 son iguales en los dos datos
hcpc <- HCPC(data_flr[,1:11], nb.clust=-1, method="average", graph = F)
table <- table(hcpc$data.clust$clust)
names(table) <- c("Cluster 1", "Cluster 2", "Cluster 3")
table## Cluster 1 Cluster 2 Cluster 3
## 39 57 28# con el analisis PCA con los 3 componentes principales
hcpc.comp <- HCPC(as.data.frame(data_flr.comp), nb.clust=-1, method="average", graph = F)
table <- table(hcpc.comp$data.clust$clust)
names(table) <- c("Cluster 1", "Cluster 2", "Cluster 3")
table## Cluster 1 Cluster 2 Cluster 3
## 39 45 40Estudiando el dendograma con los datos completos, vemos que hay dos clusters muy grandes y otro pequeño de solamente dos individuos. Debido a que son datos con pocas variables, y como hemos visto en los mapas del análisis de componentes principales no podemos decir que haya una distribución clara de los indviduos control y autistas. Como vemos en el dendograma, los individuos de la izquierda (azul y rojo) hay más concentración de individuos Control (naranja en la barra). En cambio los individuos de la derecha del dendograma (verde oscuro) hay una mayor concentración de individuos Autistas (verde claro).
dend <- as.dendrogram(hclust_data_flr_avg)
dend %>%
set("labels_col", "orange") %>% set("labels_cex", 0.6) %>%
set("labels_col", value = c("skyblue", "#FC4E07", "forestgreen"), k=3) %>%
set("branches_k_color", value = c("skyblue", "#FC4E07", "forestgreen"), k = 3) %>%
set("leaves_pch", 19) %>%
set("leaves_cex", 0.1) %>%
plot(axes=FALSE, main = "Average dendogram")
my_colors <- ifelse(data_flr$class=="Control", "orange", "green")
colored_bars(colors = my_colors, dend = dend, rowLabels = "Class")
fviz_cluster(hcpc, main = "Cluster datos completos average clust", palette = c("skyblue", "#FC4E07", "forestgreen"))
Estudiando el dendograma con los datos compilados en componentes principales, que ocurre el mismo fenómeno pero se intuye mejor separación que los datos anteriores. Ahora las ramas azules a la izquierda del dendograma tiene una mayor concentración de individuos Autistas (barra naranja). A la derecha se concentran más los individuos Control (barras naranjas).
dend <- as.dendrogram(hclust_data_flr_avg.comp)
dend %>%
set("labels_col", "orange") %>% set("labels_cex", 0.6) %>%
set("labels_col", value = c("skyblue", "#FC4E07", "forestgreen"), k=3) %>%
set("branches_k_color", value = c("skyblue", "#FC4E07", "forestgreen"), k = 3) %>%
set("leaves_pch", 19) %>%
set("leaves_cex", 0.1) %>%
plot(axes=FALSE, main = "Average dendogram")
my_colors <- ifelse(data_flr$class=="Control", "orange", "green")
colored_bars(colors = my_colors, dend = dend, rowLabels = "Class")
fviz_cluster(hcpc.comp, main = "Cluster datos completos average clust", palette = c("skyblue", "#FC4E07", "forestgreen"))
En el mapa de calor para jerarquizar tanto indviduos como variables de los datos totales, vemos que se pueden intuir en este caso dos componentes principales de variables. Las variables de la derecha se expresarian mas en autistas.
Cluster no jerárquico: K-medias
Vamos a realizar la clusterizacion de los individuos mediante la técnica de Kmeans en los datos completos y en los datos compilados. Los clusters elegidos son 3.
El resultado de la clusterización de los datos completos da una bondad del 31.6% que no es muy buena. En el caso de los datos compilados, la bondad sube hasta un 52.2%.
set.seed(12345)
options(width=3000)
km <- kmeans(data_flr[,1:11], 3, iter.max = 1000, nstart = 25)
km## K-means clustering with 3 clusters of sizes 46, 26, 52
##
## Cluster means:
## X212864_at X208901_s_at X228949_at X215372_x_at X227618_at X217142_at X207084_at X1564876_s_at X225363_at X1563327_a_at X231541_s_at
## 1 0.7020029 0.59698114 0.7220143 -0.4597533 0.7421477 -0.3083840 -0.3079812 -0.5634089 0.8015382 -0.3338598 -0.1874718
## 2 -0.8700552 -0.90187331 -0.5372643 0.9350262 -0.8058306 0.8268853 0.1338436 1.1460884 -0.8943533 0.7140844 0.9143990
## 3 -0.2505206 -0.02976545 -0.4297301 -0.1446591 -0.3111978 -0.0936955 0.0936039 -0.1231534 -0.3187707 -0.1187506 -0.3630676
##
## Clustering vector:
## [1] 1 1 1 3 1 3 3 1 1 3 3 1 3 1 1 1 3 2 3 2 1 3 1 1 3 1 1 3 1 1 1 1 1 3 1 1 1 3 3 3 3 3 1 1 1 3 3 3 2 1 1 1 1 2 1 3 1 1 1 3 2 3 2 3 3 3 2 2 2 2 1 3 1 2 3 1 2 1 2 2 3 3 2 3 2 2 3 3 3 3 1 1 3 3 1 3 3 3 1 3 1 2 3 3 3 1 3 1 2 3 2 2 2 2 3 3 3 2 2 1 3 3 3 2
##
## Within cluster sum of squares by cluster:
## [1] 326.8471 186.9409 323.5554
## (between_SS / total_SS = 31.6 %)
##
## Available components:
##
## [1] "cluster" "centers" "totss" "withinss" "tot.withinss" "betweenss" "size" "iter" "ifault"
set.seed(12345)
options(width=3000)
km.comp <- kmeans(data_flr.comp, 3, iter.max = 1000, nstart = 25)
km.comp## K-means clustering with 3 clusters of sizes 39, 40, 45
##
## Cluster means:
## Dim.1 Dim.2 Dim.3
## 1 -2.3352935 0.1129540 0.2256079
## 2 1.8539466 0.6190645 0.3783972
## 3 0.3759685 -0.6481730 -0.5318799
##
## Clustering vector:
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124
## 2 2 2 3 2 3 3 2 2 3 3 2 3 3 2 2 3 1 3 1 2 1 2 2 3 2 2 3 2 2 3 2 2 3 3 3 2 3 3 3 3 3 2 2 2 3 3 1 1 2 3 2 3 1 2 1 2 3 2 1 1 1 1 1 3 2 1 1 1 1 2 3 2 1 3 2 1 2 1 1 3 3 1 1 1 1 3 3 3 3 2 2 3 1 3 3 3 3 2 2 2 1 3 3 3 2 1 2 1 1 1 1 1 1 1 1 1 1 1 2 3 3 3 1
##
## Within cluster sum of squares by cluster:
## [1] 134.4193 133.7248 108.9038
## (between_SS / total_SS = 52.2 %)
##
## Available components:
##
## [1] "cluster" "centers" "totss" "withinss" "tot.withinss" "betweenss" "size" "iter" "ifault"
Con respecto al número de individuos presentes en cada cluster vemos que hay diferencia entre utilizar los datos completos o compilados. En cambio, los datos compilados distingue los mismos clusters que el método jerarquico, ya que si recordamos, los valores eran los mismos (39, 49, 45).
## Cluster 1 Cluster 2 Cluster 3
## 46 26 52## Cluster 1 Cluster 2 Cluster 3
## 39 40 45Si clasificamos los individuos según la class en cada conglomerado vemos que el primero es básicamente individuos Control, el segundo individuos con Autismo y el tercio una mezcla.
## [[1]]
## [1] Control Control Control Control Control Control Control Control Control Control Control Control Control Control Control Control Control Control Control Control Control Control Control Control Control Control Control Control Control Control Control Control Control Control Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism
## Levels: Autism Control
##
## [[2]]
## [1] Control Control Control Control Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism
## Levels: Autism Control
##
## [[3]]
## [1] Control Control Control Control Control Control Control Control Control Control Control Control Control Control Control Control Control Control Control Control Control Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism
## Levels: Autism ControlEn el caso de los datos compilados podemos ver el mismo patrón pero con la salvedad que en este caso el cluster 1 hay una mayoria de individuos autistas y el cluster 2 controles.
## [[1]]
## [1] Control Control Control Control Control Control Control Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism
## Levels: Autism Control
##
## [[2]]
## [1] Control Control Control Control Control Control Control Control Control Control Control Control Control Control Control Control Control Control Control Control Control Control Control Control Control Control Control Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism
## Levels: Autism Control
##
## [[3]]
## [1] Control Control Control Control Control Control Control Control Control Control Control Control Control Control Control Control Control Control Control Control Control Control Control Control Control Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism Autism
## Levels: Autism ControlAl realizar los mapas de los conglomerados, vemos que los datos totales no separan del todo bien los conglomerados.
set.seed(12345)
kmed <- pam(dist_data_flr, k =3, diss= T)
plot(kmed, which.plots = 1, main = "kmed datos totales")
En cambio en el mapa de los datos compilados se aprecia una mejor separación entre conglomerados. Además tiene menos datos negativos en el analisis Silhoutte que interfieren en clusterizacion de forma negativa.
set.seed(12345)
kmed <- pam(dist_data_flr.comp, k =3, diss= T)
plot(kmed, which.plots = 1, main = "kmed datos compilados")
Cuando revisamos los escalamientos multidimensionales, no observamos una mejora en los datos totales y en los datos compilados.
plot(data_flr.pca$ind$coord, col = km$cluster, xlab = "components 1", ylab = "components 2", pch = 19, main = "No escalado", frame.plot = F, ylim = c(-5, 5))
abline(v = 0, h = 0, lty = 2)
mds.data_flr <- cmdscale(dist_data_flr, k = 3, eig = T)
plot(mds.data_flr$points, col = km$cluster, xlab = "components 1", ylab = "components 2", pch = 19, main = "metrico", frame.plot = F, ylim = c(-5, 5))
abline(v = 0, h = 0, lty = 2)
## initial value 23.937701
## iter 5 value 20.006185
## final value 19.704767
## convergedplot(mds.data_flr$points, col =km$cluster, xlab = "components 1", ylab = "components 2", pch = 19, main = "no metrico", frame.plot = F)
abline(v = 0, h = 0, lty = 2)
Vemos que el escalado multiple no beneficia la separacion de individuos.
plot(data_flr.pca$ind$coord, col = km.comp$cluster, xlab = "components 1", ylab = "components 2", pch = 19, main = "No escalado", frame.plot = F)
abline(v = 0, h = 0, lty = 2)
mds.data_flr <- cmdscale(dist_data_flr.comp, k = 3, eig = T)
plot(mds.data_flr$points, col = km.comp$cluster, xlab = "components 1", ylab = "components 2", pch = 19, main = "metrico", frame.plot = F)
abline(v = 0, h = 0, lty = 2)
## initial value 14.912937
## iter 5 value 12.706154
## final value 12.641103
## convergedplot(mds.data_flr$points, col =km.comp$cluster, xlab = "components 1", ylab = "components 2", pch = 19, main = "no metrico", frame.plot = F)
abline(v = 0, h = 0, lty = 2)